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We discuss sensitivity to initial conditions in a model for avalanches in granular media displaying 
self-organized criticality. We show that damage, due to a small perturbation in initial conditions, 
does not spread. The damage persists in a statistically time-invariant and scale-free form. We 
argue that the origin of this behavior is the Abelian nature of the model, which generalizes our 
results to all models with Abelian properties, including the BTW model and the Manna model. 
An ensemble average of the damage leads to seemingly time dependent damage spreading. Scaling 
arguments show that this numerical result is due to the time lag before avalanches reach the initial 
perturbation. 

PACS numbers: 89.75.-k, 89.75. Da, 05.65+b, 45.70.Ht, 45.70.Vn. 



There is evidence that the dynamics of a pile of rice 
may display self-organized criticality (SOC) In care- 
ful experiments where elongated rice grains were slowly 
dropped between two glass plates, Frette et al. found 
scale- free behavior in a rice pile [2|. In a slowly driven 
pile, the angle of repose evolves to a stationary state and 
the behavior of the system is dominated by a scale-free 
avalanche size density and punctuated equilibrium. This 
punctuation causes SOC systems to be highly non-linear, 
a single grain added at one end can result in an avalanche 
propagating through the entire system. However, SOC 
models typically have stick-slip dynamics Q and it has 
yet to be established whether they allow the non-linearity 
to manifest itself in the form of sensitivity to initial con- 
ditions, as it does in chaotic systems. 

We have studied damage spreading in a simple one- 
dimensional granular model known as the Oslo model, 
which exhibits SOC It describes a number of slowly 
driven granular systems and belongs to the same uni- 
versality class as a model for interface depinning in a 
random medium and the Burridge-Knopoff train model 
for earthquakes 0, IE • The Oslo model has largely re- 
sisted efforts for an analytic solution, the few exceptions 
have been the exact enumeration of the number of recur- 
rent configurations [jj, the mapping of the model to the 
quenched Edwards- Wilkinson equation in the continuous 
limit 0, 0, and the transition matrix results 01 an d 
operator algebra recently developed for the Oslo model 
[llj. In this letter, we add to its analytical description 
by illustrating its Abelian properties |12| . 

We find that damage due to a small perturbation does 
not spread in the Oslo model as the damage is unable 
to evolve. It is possible to represent the perturbation in 
terms of commutative operators which leads to a statis- 
tically time-invariant and scale-free damage. This phe- 
nomenon is in contrast to what is seen in chaotic and 



equilibrium systems. We also address the results ob- 
tained by a previous study on an ensemble avera ge o f 
the damage, which seems to contradict our findings [13J. 
In fact, we show that these results are consistent with 
ours and that the observed behavior may be derived us- 
ing simple scaling arguments. 
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FIG. 1: (a) The Oslo model of a one-dimensional granular 
pile. Grains are added at the site x = 1 next to the vertical 
wall by letting h(l, t+1) — h(l, t) + l. A grain at site x topples 
to site x + 1 if the local slope z(x,t) — h(x,t) — h(x + l,t) 
exceeds its critical slope z c (x). When z(L,t) > z c (L), a grain 
leaves the system at the open boundary, h(L, t) — > h(L, t) — 1. 
(b) The solid grain is the additional grain in the copy at time 
t . 

The model: The Oslo model is defined on a one di- 
mensional discrete lattice with L sites at positions x = 
1, 2, ... L, see Fig. ^ On the left, the boundary is a ver- 
tical wall, and on the right, the boundary is open. The 
height of grains at site x and time t is denoted h(x, t), the 
local slope is defined as z(x, t) = h(x,t) — h(x+l,t). Each 
site has a critical slope, z c (x), which takes the values 1 
or 2 with equal probability. At each time step a single 
grain is added to the site at x = 1. If the local slope at 
any site, x, exceeds its critical slope, z(x,t) > z c {x), an 
avalanche is initiated. The site will relax and a grain will 
topple from site x to site x + 1, i.e. h{x, t) — > h{x, t) — 1 
and h(x+ 1, t) — > h(x + l,t) + 1. Each time a site relaxes 
its critical slope is redetermined, chosen randomly from 
the values 1 or 2. This toppling may cause sites x ± 1 to 
exceed their critical slopes, in which case these sites re- 
lax in turn. The avalanche will continue until the system 
reaches a stable configuration, when z(x,t) < z c (x) for 
all x. 
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consider a system which has been evolved to the critical 
state and make an exact copy at some time, t Q . Wc 
define h°(x,t) and h c (x,t) as the heights of the original 
and copy, respectively, and z°(x) and z c c {x) as the critical 
slopes of the original and copy, such that 



h c (x,t ) 



h°(x,t ) 
z° c (x) 



1 < x < L. 



(1) 



We then perturb the copy by adding a single grain to a 
site, i, such that h c (i,t ) = h°(i,t ) + 1, see Fig. ^b). 
This extra grain is not allowed to topple until it is toppled 
upon by an avalanche from above. The two systems are 
then evolved using exactly the same sequence of random 
thresholds {z c (x)} for corresponding sites in the origi- 
nal and copy. This ensures that the damage observed is 
purely due to the extra grain in the copy. Further justi- 
fication for this is to consider the mapping of the model 
to interface depinning as it is clear that the medium the 
interface moves through does not change as a result of 
the perturbation. 

We measure the damage, defined as 



H(t,L) = J2\h°(x,t) - h c (x,t)\. 



(2) 
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FIG. 2: Damage as a function of time for a single run, L — 
128, where the perturbed site is i = 16 and we have taken 
to = 0. Notice that there is no temporal evolution in the data, 
the damage is continually fluctuating between low and high 
values and frequently returns to a H = 1 configuration. The 
triangles indicate where the H = 1 configuration corresponds 
to a repaired configuration, where corresponding sites in the 
original and copy have relaxed an equal number of times. 



Figure [3 is a plot of damage versus time for a typical 
simulation. There appears to be no sense of temporal 
evolution in the data, the damage is continually fluctu- 
ating between high and low values, often returning to 
H = 1 where the original and copy only differ by a sin- 
gle grain. The triangles in Fig. [21 indicate which of these 
H = 1 configurations are repaired configurations where 
the extra grain in the copy is at the site which was orig- 
inally perturbed and corresponding sites in the original 
and copy have relaxed exactly the same number of times. 
In repaired configurations, the difference between the two 
systems is exactly equivalent to the initial perturbation 
at tQ. The occurrence of a repaired configuration corre- 
sponds to 'resetting the clock', meaning that the damage 



may not evolve in a single pair of systems, it is statisti- 
cally time invariant. 

This behavior emerges because the Oslo model has 
an Abelian nature, where the commutative operation is 
adding one unit of slope to a site and allowing the system 
to relax. If we add one unit of slope at site x and then 
at site y, it is equivalent to adding at site y and then x. 
The Abelian nature is alluded to by the fact that these 
simulations require a careful record of the sequences of 
random critical slopes. In this way the sequences of crit- 
ical slopes are no longer random noise, but an intrinsic 
property of the system where the values, although gener- 
ated randomly, are treated as given a priori. Note that 
the model is not Abelian in the strict sense of an Abelian 
group as there is no inverse operation. 

We introduce the notation C (t) to represent the stable 
configuration of a system at time t. The relevant opera- 
tors are the perturbation operator Pj, which simply adds 
one unit of slope to site i and Pi , which adds one unit of 
slope to site i and allows the system to relax if necessary. 
Thus, Pi represents a mapping within the configuration 
space of C(t). The evolution of C(t) can be expressed by 
an evolution operator T = Pi, that is 



C{t+\)=TG(t). 



(3) 



The operators Pi and T obey the commutation relation 



P l TC(t) =TPiC(t), 



(4) 



which we prove elsewhere |l4j . Hence, adding slope to 
a system and then evolving it leads to the same config- 
uration as evolving it and then adding slope, as we may 
always move the Pi operator on the right hand side of 
Eq. J3J to the left of the T operator. 

If we begin with a system at a point C(0) in configura- 
tion space, then the repeated operation of T defines the 
chain of configurations which C(t) will pass through dur- 
ing its evolution. Consider a succession of points C(t), 
which is related to C(t) by the relation (7(0) = PiC(0). 
The points C(t) and C{t) will be connected through the 
operator Pi at every step of the evolution of the two sys- 
tems, i.e. C(t) = P l C(t), for all t. 

The perturbation we studied in the simulations was 
the addition of a single grain and not slope. Adding a 
single grain to site i will change the slopes such that 
z c (i, t ) = z°(i, t ) + 1 and z c (i - 1, to) = z°(i - 1, t ) - 1. 
To analyze this situation, we start with a master system 
C m (0) and derive from this two other systems C°(0) and 
C c (0) through the relations 



C°(0) 
C c (0) 



= P 4 -iC m (0) 
= P l C m {0). 



(5) 



It is straightforward to verify that the configurations 
C°(0) and C c (0) differ by one grain at site i, thus repro- 
ducing the original and perturbed systems of our simu- 
lations with the perturbative grain placed at site i, see 
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FIG. 3: The relationship between the master, C m (0), orig- 
inal, C°(0), and copy, C c (0) configurations at time to = 0, 
see Eq. . The shaded grains are those added to the master 
configuration to produce configuration of the original system. 
The solid grain is the additional grain added to the configu- 
ration of the copy at time to. 



Fig. [3J Each of these points start a chain of stable con- 
figurations, linked by the operator T . After an avalanche 
has reached site i, each configuration in the master chain 
links to those belonging to the original and copy by the 
operators Pj_i and Pi, respectively, due to Eq. J3J. Note 
that the operator Pi does not have a unique inverse so 
there is no direct path between the original and copy. 

The important result from this is that damage cannot 
spread. Spreading implies that the small amount of dam- 
age at to leads to a little more damage at to + 1, and so 
on until the damage saturates the system, which is the 
case for a chaotic system. The Abelian nature means 
that the value of the damage at any time is independent 
of when the perturbation took place. Hence, spreading is 
not possible. However, the damage does not remain con- 
stant or decay, as in an equilibrium system, because the 
damage is exactly that due to the avalanches which would 
be caused by adding slope to different sites in the master 
system at that time. This leads to a damage size den- 
sity that is related to the avalanche size density, as easily 
recognized by considering the special case of perturbing 
site i = 1. As the avalanche size density is scale free, we 
find that the damage size density is scale free also. In an 
infinite system the damage may become arbitrarily large, 
yet it will frequently return to an H = 1 configuration! 
In other words such a system may be considered as lying 
on the edge of chaos 0, 0] . 

The damage in a single pair of systems is statistically 
time invariant, yet a previous study has found that the 
ensemble averaged data is not ^| . The ensemble average 
of damage over N runs, gives the average damage as a 
function of time 



1 N 

(H)(t,L) = -'£H j (t,L,i j ), 



(0) 



where z and (3 are exponents to be determined and G(x) 
is constant for x <C 1 and proportional to x~ z for x 3> 
1. The apparent time dependence arises from the fact 
that the perturbed site is not allowed to relax until it is 
relaxed upon from above. This forces all the systems into 
a repaired configuration at the start of the simulation and 
the average damage increases over time as more systems 
have avalanches which reach the perturbed site. 

In the case where the perturbative grain was placed 
on a random site for each system in the ensemble, scal- 
ing arguments may be used to derive the temporal evo- 
lution of (H)(t,L). The derived equation agrees well 
with the simulation data. First, we calculate how long 
avalanches take to reach the perturbed site, which we 
denote as site i. The linear avalanche size, I, is known 
to be related to the avalanche size, s, by s cx l D ±_ where 
D is the avalanche fractal dimension, D w 2.25 [5j. As- 
suming s = l D , the probability Pi(l,L)dl of having an 
avalanche with linear length in the range I — > I + dl, 
obeys Pi(l,L) dl = P s (s,L) ds, where P s (s,L)ds is the 
probability of having an avalanche of size s in the range 
s — > s + ds, given by the scaling ansatz 



P s (s,L)ds = s T Gs{-£pj ds, 



(8) 



where Q s {x) is constant for x <C 1 and decays rapidly for 
i»l, and r is the avalanche exponent, t « 1.55 |5j. We 
find 



Pl{l,L) dl = r u Qi [ - ] dl, 



(9) 



where Qi (x) is constant for i<1 and decays rapidly for 
x — ► 1. Note that we have used the scaling relation r = 
2 — 1/1? which is derived from the fact that (s) = l[T^|. 

This result allows us to calculate the probability of hav- 
ing an avalanche of linear size larger than some distance 
X, which we denote <j>(X). 



4>(X) = / Pi(l,L) dl cx X 
Jx 



l-D 



(10) 



Hence, we may expect to have an avalanche of size I > X 
within the timescale 



t 



1 



X D-1 = X X. 



(11) 



where we have used the scaling relation x = D — 1, re- 
lating the roughness exponent \ to the fractal dimension 
D 0. We use Eq. (fTTfl to obtain an ansatz for (X)(t, L), 
the average linear distance reached by the avalanches in 
a time t, 



where Hj(t, L,ij) is the damage from a single run and 
the variable ij is the site of perturbation for the jth run. 
In Ref. [13J it was found that for ij — L/2 Vj 



(H)(t,L)=t z 



(7) 



(X)(t,L)=t-g 



(12) 



where ti, is the timescale after which the avalanches can 
be expected to have spanned the entire system and g(x) 
is constant for i< 1 and proportional to x~~* for i>1 
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which ensures that (X)(t,L) < L for all t. By inserting 
X = L into Eq. we see that oc L x , which leads to 
(X)(oo, L) = L as expected. 

The average damage as a function of time and system 
size, (H)(t, L), may therefore be expressed as 

(H)(t,L) = L a F H ^{t,L) + (l-F H ^L)) 

« L a F H ^{t,L) forL>l, (13) 

where Fiijti(t, L) is the fraction of systems in an H ^ 1 
damaged configuration at time t, L a is the mean amount 
of damage of a system in the damaged configuration and 
the last approximation is true for large L. If the positions 
of the perturbed sites are distributed uniformly among 
the systems in the ensemble, we expect 

(fl)(t 1 X)=L-WM=t?G(^) > (14) 

where G{x) — x~x~ g(x) and we have taken t = 0. Scal- 
ing arguments, taking care to use the scaling relation 
t = 2 — 2 /D for a bulk driven system ! yields a = 1 . 
This is consistent with our measurement of a » 1 and 
so, putting x = 5/4 and a = 1 into Eq. Ijl4|l . we have 

(H)(t,L)=t OS0 G( J ^). (15) 




FIG. 4: By plotting t -0 ' 80 (H)(t, L) versus the rescaled time 
x — L~ 1 ' 25 t the data for the ensemble average collapses onto 
a single well defined curve G(x), see Eq. 11511 . This is shown 
for system sizes L — 64, 128, 256 and 512. The number of 
systems in each ensemble is 10, 000. 



A data collapse of the data using these values is shown 
in Fig. It is in good agreement with Eq. I|15|l . thus 
supporting our explanation for the appearance of a time 
dependence in (H)(t,L). However, this is only a crude 
calculation as there are many effects we have ignored. 
For instance, we have only calculated the average time for 
avalanches to reach the perturbed site. There is actually 
a distribution of times and this will contribute to the time 
dependence of (H)(t, L). In conclusion, we have analyzed 
damage spreading in the Oslo model, showing that dam- 
age is unable to evolve in a perturbed system. The dam- 
age is statistically time invariant and scale free and thus 
allows for arbitrarily large values in infinite systems. This 
phenomenon is due to the Abelian property of the Oslo 
model, and this generalizes our result to all other models 
with Abelian properties, including the BTW model |l(, 
the Manna model jig and the model for interface de- 
pinning in a random medium [a, uj • Thus, many of the 
classic models of SOC may be considered as lying on the 
edge of chaos [HI3. Finally, we have shown how sim- 
ulations may in fact lead to a time dependent ensemble 
averaged damage and have calculated this for the case of 
random placement of the perturbative grain. 
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